### PACKAGE PREAMBLE ###

# Create a vector of package names
packages <- c("dplyr", "ggplot2", "tidyr", "tidyverse", "data.table", "readxl", "foreign", "devtools", "fixest", "sandwich", "modelsummary", "kableExtra", "xtable")

installed_packages <- packages %in% rownames(installed.packages())
if (any(installed_packages == FALSE)) {
  install.packages(packages[!installed_packages])
}

# Packages loading
invisible(lapply(packages, library, character.only = TRUE))
rm(packages, installed_packages)

options(modelsummary_format_numeric_latex = "plain")
f <- function(x) format(x, digits = 3, nsmall = 2, scientific = FALSE)

getwd()
setwd(dir = "/Users/jmc4qg/The Lab Dropbox/Jonathan Colmer/ShotSpotter_env/Journal_submissions/ReStat/Replication Folder/")

mydata_staggered <- read.dta("Analysis Data/NIBRS_analysis_staggered.dta")
mydata_staggered_55 <- filter(mydata_staggered, Ei==55 | Ei== 0)
mydata_staggered_69 <- filter(mydata_staggered, Ei==69 | Ei== 0)
mydata_staggered_115 <- filter(mydata_staggered, Ei==115 | Ei== 0)
mydata_staggered_138 <- filter(mydata_staggered, Ei==138 | Ei== 0)
mydata_staggered_148 <- filter(mydata_staggered, Ei==148 | Ei== 0)
mydata_staggered_181 <- filter(mydata_staggered, Ei==181 | Ei== 0)
mydata_staggered_229 <- filter(mydata_staggered, Ei==229 | Ei== 0)
mydata_staggered_239 <- filter(mydata_staggered, Ei==239 | Ei== 0)
mydata_staggered_260 <- filter(mydata_staggered, Ei==260 | Ei== 0)
mydata_staggered_267 <- filter(mydata_staggered, Ei==267 | Ei== 0)

DiDiT_55 <- feols(homicide_pc~temp_MP + prec_MP | ori_sample_month+ORI[tMean]+sample_month[tMean]+ORI[prec]+sample_month[prec]  + week + dow, data=mydata_staggered_55)
summary(DiDiT_55, cluster="STATE")
DiDiT_69 <- feols(homicide_pc~temp_MP + prec_MP | ori_sample_month+ORI[tMean]+sample_month[tMean]+ORI[prec]+sample_month[prec]  + week + dow, data=mydata_staggered_69)
summary(DiDiT_69, cluster="STATE")
DiDiT_115 <- feols(homicide_pc~temp_MP + prec_MP | ori_sample_month+ORI[tMean]+sample_month[tMean]+ORI[prec]+sample_month[prec]  + week + dow, data=mydata_staggered_115)
summary(DiDiT_115, cluster="STATE")
DiDiT_138 <- feols(homicide_pc~temp_MP + prec_MP | ori_sample_month+ORI[tMean]+sample_month[tMean]+ORI[prec]+sample_month[prec]  + week + dow, data=mydata_staggered_138)
summary(DiDiT_138, cluster="STATE")
DiDiT_148 <- feols(homicide_pc~temp_MP + prec_MP | ori_sample_month+ORI[tMean]+sample_month[tMean]+ORI[prec]+sample_month[prec]  + week + dow, data=mydata_staggered_148)
summary(DiDiT_148, cluster="STATE")
DiDiT_181 <- feols(homicide_pc~temp_MP + prec_MP | ori_sample_month+ORI[tMean]+sample_month[tMean]+ORI[prec]+sample_month[prec]  + week + dow, data=mydata_staggered_181)
summary(DiDiT_181)
DiDiT_229 <- feols(homicide_pc~temp_MP + prec_MP | ori_sample_month+ORI[tMean]+sample_month[tMean]+ORI[prec]+sample_month[prec] + week + dow, data=mydata_staggered_229)
summary(DiDiT_229, cluster="STATE")
DiDiT_239 <- feols(homicide_pc~temp_MP + prec_MP | ori_sample_month+ORI[tMean]+sample_month[tMean]+ORI[prec]+sample_month[prec] + week + dow, data=mydata_staggered_239)
summary(DiDiT_239, cluster="STATE")
DiDiT_260 <- feols(homicide_pc~temp_MP + prec_MP | ori_sample_month+ORI[tMean]+sample_month[tMean]+ORI[prec]+sample_month[prec] + week + dow, data=mydata_staggered_260)
summary(DiDiT_260, cluster="STATE")
DiDiT_267 <- feols(homicide_pc~temp_MP + prec_MP | ori_sample_month+ORI[tMean]+sample_month[tMean]+ORI[prec]+sample_month[prec] + week + dow, data=mydata_staggered_267)
summary(DiDiT_267, cluster="STATE")


n_55 <- nobs(DiDiT_55)
n_69 <- nobs(DiDiT_69)
n_115 <- nobs(DiDiT_115)
n_138 <- nobs(DiDiT_138)
n_148 <- nobs(DiDiT_148)
n_181 <- nobs(DiDiT_181)
n_229 <- nobs(DiDiT_229)
n_239 <- nobs(DiDiT_239)
n_260 <- nobs(DiDiT_260)
n_267 <- nobs(DiDiT_267)

temp_MP_coef_55 <- coef(DiDiT_55)["temp_MP"]
temp_MP_coef_69 <- coef(DiDiT_69)["temp_MP"]
temp_MP_coef_115 <- coef(DiDiT_115)["temp_MP"]
temp_MP_coef_138 <- coef(DiDiT_138)["temp_MP"]
temp_MP_coef_148 <- coef(DiDiT_148)["temp_MP"]
temp_MP_coef_181 <- coef(DiDiT_181)["temp_MP"]
temp_MP_coef_229 <- coef(DiDiT_229)["temp_MP"]
temp_MP_coef_239 <- coef(DiDiT_239)["temp_MP"]
temp_MP_coef_260 <- coef(DiDiT_260)["temp_MP"]
temp_MP_coef_267 <- coef(DiDiT_267)["temp_MP"]

prec_MP_coef_55 <- coef(DiDiT_55)["prec_MP"]
prec_MP_coef_69 <- coef(DiDiT_69)["prec_MP"]
prec_MP_coef_115 <- coef(DiDiT_115)["prec_MP"]
prec_MP_coef_138 <- coef(DiDiT_138)["prec_MP"]
prec_MP_coef_148 <- coef(DiDiT_148)["prec_MP"]
prec_MP_coef_181 <- coef(DiDiT_181)["prec_MP"]
prec_MP_coef_229 <- coef(DiDiT_229)["prec_MP"]
prec_MP_coef_239 <- coef(DiDiT_239)["prec_MP"]
prec_MP_coef_260 <- coef(DiDiT_260)["prec_MP"]
prec_MP_coef_267 <- coef(DiDiT_267)["prec_MP"]

get_se <- function(model, var_name) {
  return(sqrt(vcov(model, cluster = "STATE")[var_name, var_name]))
}

temp_MP_se_55 <- get_se(DiDiT_55, "temp_MP")
temp_MP_se_69 <- get_se(DiDiT_69, "temp_MP")
temp_MP_se_115 <- get_se(DiDiT_115, "temp_MP")
temp_MP_se_138 <- get_se(DiDiT_138, "temp_MP")
temp_MP_se_148 <- get_se(DiDiT_148, "temp_MP")
temp_MP_se_181 <- get_se(DiDiT_181, "temp_MP")
temp_MP_se_229 <- get_se(DiDiT_229, "temp_MP")
temp_MP_se_239 <- get_se(DiDiT_239, "temp_MP")
temp_MP_se_260 <- get_se(DiDiT_260, "temp_MP")
temp_MP_se_267 <- get_se(DiDiT_267, "temp_MP")

prec_MP_se_55 <- get_se(DiDiT_55, "prec_MP")
prec_MP_se_69 <- get_se(DiDiT_69, "prec_MP")
prec_MP_se_115 <- get_se(DiDiT_115, "prec_MP")
prec_MP_se_138 <- get_se(DiDiT_138, "prec_MP")
prec_MP_se_148 <- get_se(DiDiT_148, "prec_MP")
prec_MP_se_181 <- get_se(DiDiT_181, "prec_MP")
prec_MP_se_229 <- get_se(DiDiT_229, "prec_MP")
prec_MP_se_239 <- get_se(DiDiT_239, "prec_MP")
prec_MP_se_260 <- get_se(DiDiT_260, "prec_MP")
prec_MP_se_267 <- get_se(DiDiT_267, "prec_MP")

weights <- c(n_55, n_69, n_115, n_138, n_148, n_181, n_229, n_239, n_260, n_267)

temp_MP_aggregated_coef <- sum(weights*c(temp_MP_coef_55, temp_MP_coef_69, temp_MP_coef_115, temp_MP_coef_138, temp_MP_coef_148, temp_MP_coef_181, temp_MP_coef_229, temp_MP_coef_239, temp_MP_coef_260, temp_MP_coef_267) / sum(weights))
temp_MP_aggregated_se <- sqrt(sum(weights*c(temp_MP_se_55, temp_MP_se_69, temp_MP_se_115, temp_MP_se_138, temp_MP_se_148, temp_MP_se_181, temp_MP_se_229, temp_MP_se_239, temp_MP_se_260, temp_MP_se_267)^2) / sum(weights))
temp_MP_aggregate_t <- temp_MP_aggregated_coef/temp_MP_aggregated_se

prec_MP_aggregated_coef <- sum(weights*c(prec_MP_coef_55, prec_MP_coef_69, prec_MP_coef_115, prec_MP_coef_138, prec_MP_coef_148, prec_MP_coef_181, prec_MP_coef_229, prec_MP_coef_239, prec_MP_coef_260, prec_MP_coef_267) / sum(weights))
prec_MP_aggregated_se <- sqrt(sum(weights*c(prec_MP_se_55, prec_MP_se_69, prec_MP_se_115, prec_MP_se_138, prec_MP_se_148, prec_MP_se_181, prec_MP_se_229, prec_MP_se_239, prec_MP_se_260, prec_MP_se_267)^2) / sum(weights))
prec_MP_aggregate_t <- prec_MP_aggregated_coef/prec_MP_aggregated_se


agg_coefs <- c(temp_MP_aggregated_coef, prec_MP_aggregated_coef)
agg_ses <- c(temp_MP_aggregated_se, prec_MP_aggregated_se)
agg_N <- 8501265

get_stars <- function(tval) {
  if (abs(tval) > 2.576) return('***')
  if (abs(tval) > 1.96) return('**')
  if (abs(tval) > 1.645) return('*')
  return('')
}

# Initialize LaTeX table string
latex_table <- "\\begin{table}\n\\centering\n\\begin{tabular}[t]{lc}\n\\toprule\n"

# Column names
latex_table <- paste0(latex_table, "  & (1) \\\\\n")

# Insert coefficients and standard errors
for(i in 1:length(agg_coefs)) {
  tval <- agg_coefs[i] / agg_ses[i]
  stars <- get_stars(tval)
  latex_table <- paste0(latex_table, sprintf("aggregated_coefficient_%d & %.6f%s \\\\\n", i, agg_coefs[i], stars))
  latex_table <- paste0(latex_table, sprintf(" & (%.6f) \\\\\n", agg_ses[i]))
}
# Insert number of observations
latex_table <- paste0(latex_table, "\\midrule\nNum.Obs. & ", agg_N, " \\\\\n")

# Close table
latex_table <- paste0(latex_table, "\\bottomrule\n\\multicolumn{2}{l}{\\rule{0pt}{1em}* p $<$ 0.1, ** p $<$ 0.05, *** p $<$ 0.01}\\\\\n\\end{tabular}\n\\end{table}")

# Write to .tex file
write(latex_table, "Figures and Tables/Table_3/column_5.tex")
